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Abstract 



O 

o 

Oh 1 In this paper a Monte Carlo generator, GaGaRes, is presented which can be 

• used to describe two-photon resonance production in e + e~ collisions. The 

program can generate the five lowest-lying C=+l meson states of any qq 
{Sj | combination together with the outgoing electron and positron. The depen- 

dence on the photon virtualities Q\ and Q\ is fully taken into account. The 
program also generates the density matrices of the resonance, which form an 
essential tool in the description of the decay of the resonances. Furthermore, 
^ . the program is applicable for all tagging conditions. 

Q\ 

o 



PROGRAM SUMMARY 

Ph! Title of program: GaGaRes 

Qh| Program obtainable from: CPC Library 

Computer for which the program is designed and others on which it is operable: Any com- 
puter that runs under UNIX and can handle quadruple precision numbers (both real and 
complex). FORTRAN 90 offers the possibility to define own data types. This should be 
used on systems where complex numbers are not available in quadruple precision. 
Operating system under which the program has been tested : UNIX 
Programming language used : FORTRAN 90 
Number of lines: 6009 (including comments) 
Keywords: Monte Carlo, two-photon, e + e~, resonance production 

Nature of physical problem: With the advent of LEP2 higher energies for two-photon 
reactions became available with high luminosities. This makes it possible to search exper- 
imentally for heavier resonances created in two-photon collisions and also to determine 
the dependence of the two-photon cross sections on the virtualities Q\ and Q\. Moreover, 
the decay distributions of the resonances can be studied. These experimental possibilities 
make it desirable to have a program, which can simulate events as expected from our 
theoretical understanding of resonance production by two photons. 

Methods of solution: A model based on the hard scattering approach is used to describe 
the production of the resonances [[[]]. For an exact description of the decay of the produced 
resonance the density matrix is required. Weyl-van-der-Waerden spinor calculations are 
used to obtain these density matrices. Events consisting of the momenta of the resonance 
and the outgoing electron and positron are generated by Monte Carlo methods and are 
distributed according to the theoretical cross section. 
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Typical running time: Depends on the requested accuracy and the generated resonance. 
On an Enterprise450, one can produce 380 1 Sq resonances/sec and 17 3 P2 resonances/sec, 
including the corresponding density matrix. Without density matrices, one can generate 
236 3 P2 resonances/sec. These numbers are for weighted events. For unweighted events 
these rates are typically an order of magnitude lower. 



1 Introduction 

With the advent of LEP2 a new region for two-photon physics was entered. Its high 
energy and its high luminosity, resulting in better statistics, allows one to look for the 
dependence of the two-photon cross sections on the 77 centre-of-mass energy W and on 
the photon virtualities Q\ and Q\. The higher statistics also gives one the opportunity 
to look for the decay distributions of the produced resonances. 

The theoretical framework on which the program is based is obtained from the 
hard scattering approach in which the meson resonance is treated as a non-relativistic 
bound state of a quark (q) and an anti-quark (q). In order to include also the production 
of lighter resonances some modifications had to be implemented. These modifications are 
discussed in Q. The two-photon widths of the resonances that are calculated using this 
model give a reasonably good agreement with experiment, also for the lighter resonances. 

The matrix elements following from this model have been implemented in the program, 
including all dependences on W and Q\. This is in contrast to most other two-photon 
MC generators 0, |J where the Q\ dependence is ad-hoc parametrized by multiplying the 
cross section by the square of a form factor F 2 (Q\ 1 Q\). The VMD model predicts the 
pole-mass form factor 

„ „ _ My My 

F(Ql ' ^ = (My+Ql) (My + Ql) ' 

In this form factor the scaling mass My is the mass of a resonance. For the lighter 
resonances one takes V — p, for the cc resonances V — J ftp whereas for the bb resonances 
V = T. Note that this form factor causes a very strong decrease of the cross section with 
high Qfs. In it is discussed that the model used in GaGaRes results in a much weaker 
decrease of the cross section with high Qf. 

The total angular momentum of the resonance is denoted by J, its spin by S{= 0, 1) 
and its orbital angular momentum by L(= 0, 1, 2, . . . or S, P,D,...), resulting in the well- 
known spectroscopic notation ^ 25+1 ^Lj, with parity P = (— 1) L+1 and charge conjugation 
C = (-1) L+S . 

In two-photon reactions one can, in lowest order, only produce C-even resonances. 
The present version of GaGaRes describes the production of 1 S (CT), 3 Pj (J = 0, 1,2) 
(J + ) and l D 2 (2~) resonances. The resonance can be composed of light (u, d, s) quarks 
in the states 7 = 0,1, where I is the isospin, or of heavy c and b quarks. 

The outline of the paper is as follows. We start by introducing the kinematics and 
notation and present the different calculational methods to obtain the expressions for the 
matrix elements squared. Then the discussion focusses on the density matrices, which 
are required to describe the decay of the resonances. This is followed by a description of 
the two Monte Carlo generation schemes available in GaGaRes for the integration of the 
differential cross section. Finally, the structure of the program is presented, together with 
some important variables, common blocks and procedures. 
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2 Notation, Kinematics and Calculations of the Ma- 
trix Elements Squared 

The process of resonance production in two-photon physics is described by the following 
reaction 

e + (pi) + e-{p 2 ) -> e + (p 3 ) + e"(p 4 ) + Y(h) + J*{k 2 ) 
^e+(p 3 )+e-(p 4 ) + R(P). 

In this reaction the 7*'s are off-shell photons which collide to form the resonance. 
Four-momentum conservation in the vertices gives 

h=pi- p 3 , 

k 2 =p 2 - Pi, (3) 
P = fci + k 2 . 

The following parametrization of the external four-momenta in the lab-frame has been 
chosen. 



(4) 



In these formulae E b is the beam-energy. P b is the momentum of the incoming leptons. 
The polar angles #3 and 64 are defined as the angle between the incoming lepton and the 
corresponding outgoing lepton. 

It is convenient to introduce the set of invariants 



Pi = 


(E b ,0,0,-P b ), 






P2 = 


(E b ,0,0,Pb), 






P3 = 


(E 3 , \p 3 \ sin ^3 cos 3 


, \p 3 1 sin 3 sin 03, 


-|p 3 | COS0 3 ) , 


Pi = 


(E4, \pi sin 4 cos 04 


, p*4 sin 4 sin 04, 


p 4 COS 4 ) , 


P = 


(E R , \p R \ sm9 R cos0 


< R , \p R \ sin 9 R sin, 


i>R, \pr\ cos 9 r 



s 


= (P1+P2) 2 , 




s' 


= (P3+Pi) 2 , 




t 


= k 2 = -Q\ = (pj 


-Ps) 2 


f 


= k\ = -Ql = (p 2 


-Pi) 2 


u 


= (pi - Pi) 2 , 




u' 


= (P2-P3) 2 , 




w 2 


= P 2 = (h + k 2 ) 2 , 





(5) 



where Ql and Q\ are called the virtualities of the intermediate photons. These invariants 
satisfy the relation 

s + s' + t + t' + u + u = W 2 + 8m 2 e . (6) 
Using the methods of expressions for the matrix elements of reaction (Q) can be 
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(7) 



obtained^ 

M^So; Ai, A 2 , A 3 , A 4 ) = [fei, k 2 ,ji{Xi, A 3 ), j 2 (A 2 , A 4 )] , 
.M( 3 P ; Ai, A 2 , As, A 4 ) = § ([MA!, A 3 ) • j 2 (A 2 , A 4 )A; 1 • fc 2 

-Ji(Ai, A 3 ) • fc 2 j 2 (A 2 , A 4 ) • h] (W 2 + h ■ k 2 ) 

-Ji(Ai,A 3 )-j 2 (A 2 ,A 4 )tt') ) 
M( 3 Pi, A l5 A 2 , A 3 , A 4 , A B ) = § its [e*(A*),.7i(Ai, A 3 ), j 2 (A 2 , A 4 ), fc 2 ] 

+fe [e*(A fl ), j 2 (A 2 , A 4 ),ii(A 1 , A 3 ), h}) , 
M( 3 P 2 ; A 1; A 2 , A 3 , A 4 , A fl ) =§(h- k 2 j lfl (\i, \ 3 )j 2v (\ 2 , A 4 ) 

+/cii»A^2i/ji(Ai, A 3 ) ■ j 2 (A 2 , A 4 ) 

-ki^,j 2v {\ 2 , A 4 )ji(Ai, A 3 ) • k 2 

-^(Ai, A 3 )j 2 (A 2 , A 4 ) ■ fcx) £*^(A B ), 
M( 1 J D 2 ; Ai, A 2 , A 3 , A 4 , A fl ) = §e*^(X R )k llx k 2v 

£ [fci,fc2,ji(Ai, A3), j' 2 (A 2 , A 4 )] . 

In these formulae Aj represents the helicity of the associated external particle. The lepton 
currents j\ and j 2 are defined by 

Ji (Ai, A 3 ) = v Xl (j)i)Yvx 3 (p3), 
J 2 (A 2 , A 4 ) = ^4(^)7^2(^2), 

satisfying current-conservation: 

ji(A m , A n ) • ^ = (i = 1,2). (9) 

The polarization vector of the outgoing spin-1 resonance with helicity A B is £^*(A B ) 
whereas e^ u *(\r) is the polarization tensor of the outgoing spin-2 resonances with he- 
licity \r. We use the polarization vectors as defined in (3.15) of f|. With this choice, the 
spin-2 polarization tensors can be constructed from the spin-1 polarization vectors using 
the Clebsch-Gordan coefficients 

e*^(±2) = e*"(±l)e w, (±l), 

e*^{±l) = ±^ (e*"(±l)e* v (0) + e^(0)e**'(±l)) , (10) 
£ *^(0) = (-e*^(+l)e* u (-l) + 2e*»(0)e* v (0) - e*»{-l)e* v (l)) . 

Finally, the constants c\, . . . , C5 are given by 

ci = g , c 2 = 4 9l /W, c 3 = 2VQgi, c 4 = aV3W 9i , c 5 = 8V30g 2 , (11) 

where we have introduced 




16e 2 e^>(0)|c* 
9% (s + s' + u + u'-SmD^^lW [ ' 

The values used for the fractional charges e 2 , and for (the derivatives of) the radial part 
of the wave functions in the origin, |i?(0)|, can be found in [|TJ. 

The above formulae lead to expressions and numerical results for the total matrix 
element squared and the density matrices. We thus follow a calculational strategy that 

1 For the Minkowski metric we use (1,-1,-1,-1). The totally antisymmetric Levi-Civita tensor is 
denned by £0123 = +1. Also the notation e\p,q,r,s] = e^ l/Xp p tl q v r\s p is introduced. 



differs from the usual one ||, in which the expression for the total cross section is split 
into three parts: two density matrices for the production of the virtual photons from 
the external leptons and a cross section for the process 7*7* — > R, which in turn can be 
written in the well-known form 



+2p 1 °pi + asT + pfpfass + 2|p+-p 2 + -|r rr cos(20) (13) 

R\n+°n+°\T d 3 $ 3 d 3 ff 4 
-o\Pi P2 \ T Tscos{(p)j -g^-ET, 

where <fi is the angle between the two scattering planes in the 77 rest frame. 

In GaGaRes there are three methods implemented for the calculation of the square of 
the matrix element summed over the helicities of the external particles. 

In the first method we use an expression for the total matrix element squared in 
terms of the invariants introduced in (|5|). These expressions have been obtained using 
FORM 0. In the calculation we have introduced the tensors formed by the product of a 
lepton-current with its complex conjugated, summed over all helicities 



£f = EWftAi, A 3 )jT(Ai, A 3 ) = + + \tgT), 
£f = Ea 2 ,a 4 J 2 m (A 2 , A 4 )jT(A 2 , a 4 ) = mpl+pip v 2 + 5<V)- 



(14) 



We also use the standard summation rules for the polarization vectors/tensors of massive 
particles 

Ex R e flu (X R )e*^(X R ) = l(V tip V ua + V^V up )-lV flu V pa . 1 ] 

The results of this exercise are collected in appendix [A|. The numerical calculation of the 
total cross section involving these expressions is the fastest method in GaGaRes. 

We now turn to the two other calculational methods, which are options in the program. 
They are based on the evaluation of helicity amplitudes in the Weyl-van-der-Waerden 
(WvdW) spinor formalism 0, [5]]. These methods are slower than the method using the 
expression in terms of invariants, but they allow a straightforward construction of density 
matrices for the produced resonances. In the WvdW calculations we use the notation and 
conventions introduced in ||]. Furthermore the Levi-Civita tensor can be written in the 
WvdW formalism as 

£ AWBXCYDZ = 4l \^AB £ CD £ WY^XZ ~ £ AC £ BD^WX^YZ 

where Eab is the metric tensor for the Weyl spinors, given by 

^ = "ab = e AB = e AB = (17) 

This metric tensor also defines the (antisymmetric) inner product of two Weyl spinors 

< pk >= e AB p A k B = PAk A = - < kp>, <pp >= 0. (18) 

In the WvdW formalism each Minkowski four-vector k^ can be related to a WvdW 
bispinor 

k> = {k°, k\ k 2 , k 3 ) - K AB = = + £ £ + ^\ (1Q) 
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As a consequence the inner products of two Minkowski four-vectors can be written as 

k-p= l -K AB P AB = l -K AB P ABl k 2 5 c B = K AB K AC . (20) 

A WvdW bispinor corresponding to a time-like (k 2 > 0) or light-like (k 2 = 0) Minkowski 
four-vector k^ 1 = (k°, \k\ sin#cos0, \k\ sin#sin0, \k\ cos8) can be decomposed into two 
Weyl spinors 

K AB = E K i,A K i,B, (21) 
i=l,2 

with 

„,^v^( e ";: o ? sf ), ^=\ra(_ e ti f ). (22) 

From these formulae one sees that for light-like four-vectors the WvdW expressions sim- 
plify dramatically as for these four-vectors the Weyl spinors with subindex 2 vanish. For 
these four-vectors the corresponding WvdW bispinor can be written in a simple dyad 
form. 

The usual Dirac spinors for helicity ±| particles and antiparticles can also be expressed 
in the above Weyl spinors 

U + = (— K A , K l ^), M_ = (« x , K^), = (ft^, — K 2 ^), W_ = (/Cg , (24) 

Using 

,-(4^), (25) 
the currents (S) can also be written in the WvdW formalism, for positrons one has 



Jf 5 (++) 


= 2 


(P3 


)f(Pl 


)f + (Pi 


)f(P3 


)fj 


J? s (+-) 


= 2 


(P3 


)f(pi 


)f " (Pi 


)?(P3 


)f 


J? s (-+) 


= 2 


_(P3 


)f(pi 


)f - (Pi 


)f(P3 


)f 




= 2 


(Pi 


)f(P3 


)f + (P3 


)f(Pl 


)f 



(26) 



The currents for the electron line can be obtained from these currents by making the 
substitutions p\ <-> p 2 and p 3 p 4 . 

It can be checked that these WvdW bispinors satisfy current-conservation, which reads 
in WvdW notation 

J AB (X m ,Xn)K hAB = (i = l,2). (27) 

Finally, the polarization tensors for the massive spin-2 resonances can be constructed 
from the spin-1 polarization vectors using the Clebsch-Gordan coefficients given in fllPf) . 
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The expressions for the amplitudes in ([7|) read in the WvdW notation 
-M(%)(Ai, A 2 , As, A 4 ) = & { J^CAi, M)K? Y K 2>BZ J$ Z {\ 2) A 4 ) 



-M( 3 P )(Ai,A 2 ,A 3 ,A 4 ) 
-M( 3 -Pi)(Ai, A 2 , A 3 , A 4 , Xr) 



eg 

lft 



^( 3 P 2 )(Ai,A 2 ,A 3 ,A 4 ,A R ) 



M( 1 J D 2 )(A 1 ,A 2 ,A 3 ,A 4 ,A fl ) 



^3)K^ X K 2 J j X J^ Y (X 2 , A 4 )| , 
7 |2i7 Jx,ab(Ai) A3) J 2 ^ jB (A 2 , A 4 ) 
— ^^i,ab(Ai) A3) J 2 (j D (A 2 , A 4 )K 2 ^ B -ft'f D | , 
^ {^2,cy(A2, \±)e*^ Y (\ R )J l BZ (\ 1 , \z)K% z 
~^2,cy(^2, X^)e* x (Xr) J 1 ; x)x(Al; A3)-ftT 2 Dy ] 
+t / [/ 1) cv(Ai, A 3 )^ y (A*) J 2iBZ (A 2 , A 4 )i^P 
— ^i,cy(^1' A3)£* cx (Ar) J 2 £, x (A 2 , A 4 )fTf yi 
^7 {2P J l y 4 B (Ai, A 3 ) J 2 ( j, D (A 2 , A 4 ) 

+-^l i AB-^2,CD^l,BF(Al, A 3 ) J 2 BF (A 2 , A 4 ) 

~K 1 ^ B J 2 ^ D (X 2 , A 4 ) J x g F (Ai, As)^^ 
"^AbA.cdCAi, A3)^ 2j sf(A2, A 4 )i^f F | 
^ D (A fi ), 

>f( 1 5 , o)(Ai,A 2 ,A 3 ,A 4 ), 



(2? 



where 



(29) 



F —ki- k 2 , 
G =W 2 + F, 
H = FG — til. 

The second and third calculational methods in GaGaRes both use the WvdW calcu- 
lation. In the second method the matrix elements are finally expressed in terms of the 
standard spinor inner products as defined in equation (^). For the derivation of these 
expressions FORM has been used. 

The third method expresses the matrix element in traces of 2 x 2 matrices. Let us 
introduce a down-matrix 



Ki = K AB , 



(30) 

and an up-matrix lO 

K^ = (K^) T , K^=K AB . (31) 

With these matrices a product of an even number of WvdW bispinors can be rewritten 
as the trace of a product of the previously defined down- and up-matrices, e.g. 



p ab k 



PabQ 



AB 

CB p cAD 
rL(j D D 



Tr 
Tr 



P^Ki 



2p ■ k, 



P1QIR1& 



(32) 



One also needs the expressions for the spin-2 polarization vectors in this method 



Tr 
Tr 

Tr 



Ph* n (±l)K^ 
PV^(O)^ 



Tr 



PV[(±1)]1 Tr \KhK±l)] 



$| (Tr [PT £I(± 1 



+Tr 
1 

VE 
+2Tr 

-Tr 



P T eJ(0) 
■Tr ' 



Tr 



Tr 



AT T eT(0) 



P^K+l] 



Tr 



(33) 



P T eJ(0)J Tr [^^(0) 
P T eJ(-l)] Tr [K^K+l) 
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As a result the matrix elements can be written as 
Mi 1 S ) 

M( 3 P ) 

M(?Pl) 

M( 3 P 2 ) 



M( 1 D 2 ) 

It turns out that this is the fastest method using the WvdW formalism, thanks to a fast 
matrix multiplication in FORTRAN90. The other method using the WvdW spinorial 
inner products merely serves as a cross-check. 

For the evaluation of the matrix element squared or of the density matrices, sum- 
mations over the helicities of external particles have to be made. However, the matrix 
element need not be evaluated for every helicity configuration. From the choice of the 
polarization vectors in |4j] it can be shown that the complex conjugate of a matrix element 
is related to the matrix element with flipped helicities: 

Spin-0 resonances : M({—Xi}) = (Hi K) M*({\i}), 

Spin-1 resonance : M({-X t }, -X R ) = (FT.; K) M*({X t }, X R ), (35) 
Spin-2 resonances : M({-X l }, -X R ) = (Eli A*) {-l) XR M*({Xi}, X r ). 

Here again the Aj's are referring to the incoming and outgoing leptons and X R is the 
helicity of the resonance. 

It should be noted that in the limit m e —>■ 0, where the WvdW calculations simplify 
dramatically, we have analytically checked that the WvdW calculation gives the same 
result for the total matrix element squared as the massless parts of the cross sections in 
appendix |X[ 

3 Density Matrices 

When one considers the decay of a resonance into some final state X, 

e + e~ -> e + e"R -> e + e~X, (36) 

one can write for the amplitude 

M = J2aP,W)A XR V XR . (37) 



_£LTr 
4tt> ±L 



J u (A 1 ,A 3 ){^ 2i J 2 T (A 2 ,A 4 ) 
-4(X 2 ,X 4 )K 2l Kl}], 
^ 7 (2^Tr[j u (A 1 ,A 3 )J 2 T (A 2 ,A 4 )" 
-GTr [j u (A 1; X 3 )Kl] Tr [j 2i (A 2 , A 4 )^J 
^ (tTr [J 2 T (A 2 , A 4 ) {£j(A K ) Ji T (Ai, X 3 )K 2l 
-K 2l Jl(X 1 ,X3)el(X R )} 
+t'Tr [f 1 (X 1 , A 3 ) {el(X R )4(X 2 , X,)K n 
-K u f 2 (X 2 ,X A )el(X R )}]), 

(2FTr [ji(Ai, X 3 )e* [i (X R )J 2 (X 2 , A 4 ) 



8tt- 

+Tr 
-Tr 
-Tr 



lc 



■Tr 



K 1 r £^(A R )J 2 T (A 2 ,A 4 ; 
K 2 T £| i (A i? )J 1 T (A 1 ,A 3 

n 



•Aj.(Ai, A 3 ) J 2 (A 2 , A 4 ) 
Tr [j u (Ai,A 3 )A: 2 T 
Tt\J 21 (X 2 ,X 4 )kI 



KleUX^K^Mi'So). 



(34) 
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A\ R describes the (two-photon) production of a resonance with helicity A#. This is 
our previously discussed matrix element. T>x R describes the decay of a resonance with 
helicity Xr into the final state X. The quantity £(P, W) represents the propagator of the 
resonance and numerical factors. The total matrix element still depends on the external 
four- momenta. No implicit integrations over outgoing momenta have been carried out. 
The square of the total matrix element is then given by 

Y.\M? = UP,M)\* £ A XR yV XR y R = \aP,M)fTr(AV*). (38) 

x R' X 'r 

Here is the summation over the initial and final states and the helicities of the in- 
termediate state. The quantities A\ R y R and T>\ R y R are A\ R A* x i and V\ R V* X , , summed 
over the helicities of all particles but the resonance. They are the density matrices for 
the production and decay of the resonance. For a spin J resonance the density matrix is 
formed by a (2 J + 1) x (2 J + l)-matrix. 

GaGaRes calculates for every generated event the normalized density matrix for the 
production of the resonance 

n _ T.Xj A(\j, \ R )A*(\j, X' R ) 

where the summation in the numerator is over the electron and positron helicities and in 
the denominator the summation is over all helicities. 

The density matrix for the production only depends on the production process and 
can therefore be calculated in GaGaRes. If one knows the mechanism for the decay of 
the resonance into a state X, also that decay density matrix could be incorporated in 
GaGaRes. This is beyond the scope of the present investigation. 

From the definition of p one can immediately derive two properties of this density 
matrix: 

• p is hermit ian : p — p'; 

• The trace of the density matrix equals 1 : Tr (p) = 1. 

However, from the relation between complex conjugated matrix elements and spin 
flipped matrix elements (|35| ) and the definition for the normalized density matrix (|39|). 
one can also derive an additional symmetry. For spin-1 resonances it reads 

P\ R \' R = P-\' R -\ R - (40) 

Px R x> R = (-1) Xr+X '« P -x> r -x r . (41) 

This symmetry, the hermiticity and the fact that the trace is one, reduce strongly the 
number of independent elements in the density matrix. For a spin-1 resonance the density 
matrix has one real diagonal element and two complex off-diagonal elements, e.g. P-i,i 
and p-1,0- For a spin-2 resonance one has two real diagonal elements and 6 independent 
complex off-diagonal elements. The fact that the diagonal elements are real follows from 
hermiticity. 

In GaGaRes only the relations in (|35|) have been used to reduce the amount of cal- 
culations. However, checks have been performed that the density matrices are hermitian 
and that equations (|40D and ([|l]) hold. 
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The density matrices are computed for every event in the resonance rest frame (RRF), 
as this is the most natural scheme for the description of the resonance decay. The 
positive z-axis, the quantization axis, is defined by the boost direction. The Lorentz- 
transformation from the lab-frame to this RRF is performed by subsequently an azimuthal 
rotation, a polar rotation and a boost. The rotation brings the momentum of the reso- 
nance in the lab frame along the z-axis and the boost then brings the resonance to rest. 
From an experimental point of view this choice of polarization axis is also favoured as it 
is straightforward to reconstruct the boost direction from the detected final state X. 

Some results on the calculation of density matrices using GaGaRes and a discussion 
thereof can be found in || . 



4 MC Generation Schemes 

The parametrization of the external four-momenta in the lab-frame is given in (||). From 
now on we will use the following shorthand notations 

Si = sin is _ 
d = cos 9i. 

Additionally we introduce 6*34, the angle between the two outgoing leptons, given by 

C34 = COS #34 = -C3C4 - S3S4 COS(0 3 - 04 + 71"). (43) 

The general expression for the cross section of a 2 — > 3 process reads 

a = f (2tt)^ (p! + p 2 - p 3 - p A - P) 1 J E \M\ 2 

d% d% d 3 P 1 1 

(2tt) 3 2_E 3 {2tt) 3 2E 4 (2tt) 3 2E r - 

Performing the integration over P and rewriting the remaining delta function gives 



a 



J ( ] 8sE 3 eJi-^ (45) 



-2|p 3 ||p4| cos 034 - M 2 )d 3 p 3 d 3 p4- 
Solving the delta function for E 3 yields 



-C{E4 - 2E b ) ± D 
2{(E A -2E b y-\p^c 2 u } : 



^3 = 0771^9 1^12^2 ( 46 ) 



where 

i =AE 2 h -AE b E A + 2ml-M 2 } 
D =|p 4 | C3 4|A|, (47) 
A 2 = e~ 4(£ 4 - 2E h fm 2 e + 4|p 4 | 2 c| 4 m2. 

The sign ambiguity in ( pf) can be solved by dividing the phase-space into two parts. The 
cross sections of the two solutions are denoted by er±, corresponding to the sign in (|46|) . 
For the two solutions we find 
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— solution 



E 
E 



4,min 



4,max 



£ m e M M 2 



2E b 



m e < E4 < E' A then — 1 < C34 < +1, 
E^ 5; E A ^ E Amax then — 1 < C34 < 



(46 



+ solution 



where 



and 



E' 



E^^min — E^ 
E/^.rnax Efo ^ ^ 
— 1 < C34 < — C a , 

AE? - Am e E b + 2m 



m e M _ M 2 



M 2 



AE b - 2m e 



C a 



\ 



A(2E b - E 4 ) 2 m 2 e - (AE 2 - AE A E b + 2m 2 e — M 2 y 
A\pl\ 2 m 2 



(49) 



(50) 



(51) 



Making a change of variables and performing the integration over E% leads to the total 
cross section 



\p3\W4 



256vr 5 s A _ 4m| | 4 £ 6 -2E A + 2 



\P4\E 3 
\Pz\ 



-d£ 4 d 2 ft 3 d 2 ft 4 . 



C34I 



(52) 



From now on all all variables are normalized to the beam energy E b . After this 
normalization the cross section reads 



<7± 



\Ps\\P4 



256n 5 s 



1-^14-2^ + 2^,, 



■d^d^d 2 ^- 



(53) 



GaGaRes offers the user the choice between two schemes in which different sets of 
variables are generated to perform the MC integration. Scheme I, based on JTof, is the 
most efficient scheme for no-tag events, as the peaking behaviour of the matrix element 
is better described. Scheme II is more efficient for single- or double-tag events, as both 
polar angles are generated variables. 



Scheme I is described in [ 1 1 for the process 



e + e — > e + e l + l (/ = e,/i, r), 

where the outgoing electron and positron are not tagged. As we have here one particle 
less in the final state, some simplifications can be introduced. This generation scheme 
being very similar to the scheme in |10[, will not be discussed in this paper. 

The second scheme is more efficient not only in the generation of single- and double- 
tagged events, but also for the generation of events with cuts on Q 2 , where an automatic 
cut on the angles is imposed. 

The procedure is to start with the normalized expression given in (|S3|). The first 
approximation one makes is to replace the expression for the matrix element 



(54) 
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The corresponding weight factor is given by 

W, = (55) 

In this approximation it is better to put a factor tt' rather than a factor t 2 t' 2 in the 
denominator: from the formulae in appendix |A| it follows that in the limit where t and t! 
are of order m 2 , the numerator also is of order m 2 , for it can be shown that in this limit 
the invariant 02- is also proportional to m 2 . This reduces effectively the behaviour of the 
denominator to tt' . This approximation also gives the well-known logarithmic behaviour 
of the total cross section with respect to the total energy. 
The cross section now reads 

*± = \ 1 1 i / d£ 4 d 2 fi 3 d 2 fi 4 . (56) 

4tt 5 s ^Ji-^J tt' |4 - 2E 4 + 2^c 34 | 

Before performing the next approximations, let us consider the photon propagators 
more closely. The square of the invariant mass of the photon, radiated by the incoming 
positron, can be written as 

t = -2^3(1-03 + ^3), (57) 

with 

(l-\p 3 \P b \ m 2 

Now 63 m in is defined as 

A -A I ^e(l-^3) 2 ,, Q , 
3 ,mm = 3 | C3= i - ^2 ' ^ ) 

which reaches its minimum at E% = E^ max = 1 — |m e M — \M 2 and is denoted by 

I m 2 (m e M+|M 2 ) 2 

- Os.minl^^^ g(1 _ | mgM _ i M2)2 - lOUj 

For i' we find a similar expression with the subscript 3 replaced by 4. 
The second approximation now consists of the substitution 



1 



tt' 4£ 3 £ 4 (l-c 3 + e)(l-C4 + e)' 
The corresponding weight factor reads 



(61) 



= 4E 3 £ 4 (l-c 3 + g )(l-c 4 + g ) 

tt' ' y ' 

The third approximation consists of the approximation in the denominator 

I|M| 2 _ Bi+ K^ CmH1 _ B , (63 ) 

^3^4 |P 3 | 

This approximation works rather well, as for most generated events the outgoing leptons 
go (almost) back to back (c 34 ~ — 1). 

= f^- 1 = l\ E ■ (64) 

E 3 E 4 \ 2 -E, + ^c u \ 1 ' 
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After these approximations the cross section reads 



I I / d£ 4 d 2 fi 3 d 2 fi 4 



32tT 5 .< 



1 



Ami 



;i-C 3 + £)(l-C4 + £)(!- £4)' 



(65) 



First, the trivial azimuthal integrations of <p 3 and 4 over [<p m i n , 4> m ax] can be per- 
formed. The integrations over the polar angles and over E4 can be performed by using 
the standard integral 



dx 



a — x 



In 



0/ 



The results of the integrations are put into the last two weights, 



Wa =ln 



W 5t± = In 



l+e-C3, r , 



1-E 4 

.111/11 



In 



l+e— C4 . 



1+e— C4 



J 3,max < r3,min)yrA,max 



'4,min ) i 



(67) 



1 ^4,max 

while the "total cross section" is given by 

] 

Otot 



4ml 



32^ 0- 

The real total cross section can now be determined using 
o = a tot (< W 1 W 2 W 3 W i (W 5+ + W 5 



> 



(68) 



(69) 



where < WiW2WsW4(W5 + + W5-) > is defined as the average of the total weight. During 
the generation of events there is a probability W5- / (W5++W5-) that an event is generated 
in the phase space of the — solution. Especially for the heavier resonances almost all events 
will be generated in the — part of the solution. 

It is useful to give some comments on the generation of the integration variables. The 
variables are now generated by importance sampling in the following order: <p 3 , 4 , cos #3, 
cos 9a. Let 1Z be a random variable uniformly distributed over the interval [0, 1]. Now the 
azimuthal angles ^3 and $4 are generated according to 



^3 — W3,max 
(>4 = [4>4max 



J 3,min )^ + 03 

)n + <j)4 



(70) 



Variables that are distributed like the integrand of (pq) should be generated according to 



X (X (fl X ma x) 



Ch t)C ft 



K 



(71) 



where it should be noted that in the program 1 — x rather than x is generated in order 
to guarantee the best numerical stability. At this point cos #34 can be evaluated. Then 
E4 is generated. When E4 < E' A one has only a contribution from the — solution of the 
delta function. When E4 > E' A one has to check whether C34 is in the allowed interval. If 
so, one has contributions from both solutions of the delta function. 

Cuts on the generated variables are incorporated in the generation mechanism. They 
do not lead to a drop in efficiency (exactly the reason why scheme II is favoured for single- 
and double-tag reactions). For scheme II the program also offers the user the possibility 
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to impose cuts on E 3 , on the photon virtualities Qf and on c 34 . The cuts are applied by 
the following replacement 

\M\ 2 ^\M\ 2 \{6 } (72) 

where n & is a product of step functions representing the cuts. 

These cuts lead to a drop in generation efficiency. However, for the cuts on the photon 
virtualities Qf the angles of the corresponding polar angles are adapted in order to decrease 
this drop in efficiency. For a cut on the virtualities we use the following approximation 
formulae 



@i,min — 2arcsin 4 



Qi,min + ^ m e /yo\ 
\ ^ J -<b J -<3,max 



E Qi,min + 2m e (n) 

AE b sin 2 (^pj 

where ft is the minimum polar angle of the outgoing lepton corresponding to the photon 
virtuality Qi. On this photon virtuality the cut Qi > Qi t7n in is imposed. 

The average normalized density matrix is calculated by taking the weighted sum of 
the normalized density matrix and by dividing this sum by the total weight. This means 
that for every generated event in the lab system the momenta in the RRF are obtained 
which are then used to evaluate the RRF density matrix. 

The program generates at the same time unweighted events, i.e. events with a weight 
factor equal to 1. This is done by applying a hit or miss algorithm. 



5 Structure of the Program 

A complete overview of all subroutines can be found in the code-listing of GaGaRes. 

In the main program the file gamgaminput is read out. This file contains some input 
parameters. The variables that are read from gamgaminput are marked with an * in the 
following discussion. 



5.1 Modules 

A list of the most important modules in GaGaRes 

MODULE global_data 

In this module general variables are stored. Also some cuts and switches are contained 
in this module. 



Events* Number of events to be generated. 

StartNumber* Number of events to be generated to determine the 
maximum weight. 

IAcc Number of accepted events (after hit or miss integration). 

ICut Number of events after applied cuts. 

Eb* Beam energy. 

EWE* The expected maximum weight. 



Th3Min*,Th3Max* The minimum and maximum allowed polar angle of 
the scattered positron. 
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Th4Min* ,Th4Max* The minimum and maximum allowed polar angle of 

the scattered electron. 
Ph3Min* ,Ph3Max* The minimum and maximum allowed azimuthal angle of 

the scattered positron. 
Ph4Min*,Ph4Max* The minimum and maximum allowed azimuthal angle of 

the scattered electron. 
E3Min*,E3Max* The minimum and maximum allowed energy of the 

scattered positron. 
E4Min* ,E4Max* The minimum and maximum allowed energy of the 

scattered electron. 
OutputFile* The name of the file in which the generation 

information is written. 
ICalc* Determines which calculational method to use: 

= 1 : Use the matrix element in terms of invariants (I); 

=2 : Use the WvdW formalism in terms of 
inner products (II); 

=3 : Use the WvdW formalism in terms of traces of 
matrix products(III); 

=4 : Use I and II; 

=5 : Use I and III; 

=6 : Use II and III; 

=7 : Use all methods. 
IDens* Switch variable (0,1). When set to 1 the density matrix 

is calculated. 

ICut* Switch variable (0,1). When set to 1 checks for cuts on 

non-generated variables are made. Note that this switch 

variable is not automatically set. 
IBoost* Switch variable (0,1). When set to 1 the resonance 

momentum in its rest frame is set exactly to (M, 0, 0, 0) 

and <pR and Or in this frame are set to 0. 
It(p)Cut* Switch variable (0,1). When set to 1 the cut on 

t (t 1 ) is taken into account and the minimum allowed 

scattering angle is adapted. 
t(p)Min*,t(p)Max* The minimum and maximum allowed values of 

t and t'. 

IBreitWigner* Switch variable (0,1). When set to the invariant mass 

of the 77-system, W, is set to M, the mass of the generated 

resonance. When set to 1, W is 

distributed according to a Breit-Wigner shape 



i r 

7T <w-M) 2 +r 2 ' 



where T is the total decay width of the resonance. 
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As the photons have a Bremsstrahlung-like character, 
a low W is preferred (high weight). In order to 
keep the possibility to generate unweighted events, 
one can choose to generate the invariant 
masses only in the interval [Mr — nT, Mr + nT]. 

MODULE PhysCon 

Some physical constants are stored in PhysCon. Between brackets their values are 
given (All masses and energies in GeV): 

me m e (= 0.00051099906 GeV); 

Alpha a e . m . ' 1 " 



V 137.036^' 

Pi 7i (= 3.14159); 

BarnConv A conversion factor to go from GeV -2 to pb 
(= 3.89385 x 10 5 pbGeV 2 ). 
MODULE Resonance 

Contains some properties of the resonance. 
IWav* Determines the wave function of the resonance: 

= 1 : %; =2 : 3 P ; =3 : 3 P l5 =4 : 3 P 2 ; =5 : l D 2 ; 
IRes* Determines the composition of the resonance: 

= 1 : I = 0; =2 : I = 1; =3 : I = 1'; =4 : cc; =5 : bb; 
RName Contains the name of the resonance; 
RMass Contains the mass of the resonance; 
RWf 2 Contains the square (of the derivative) of the radial part 

of the wave function in the origin(see IJ); 
eq2 Contains the value for e 2 (see [|l|); 

GGWidth Contains the value of T 77 , the two-photon width; 
Width Contains the value of T, the total decay width of the resonance. 

5.2 COMMON blocks 

Some important COMMON blocks: 
COMMON FourVecs 

Contains the four-momenta of the generated particles in the lab-frame. 



PV1 


Pi- 


four-momentum of incoming 


positron; 


PV2 


P2- 


four-momentum of incoming 


electron; 


PV3 


P3- 


four-momentum of outgoing 


positron; 


PV4 


Pa- 


four-momentum of outgoing 


electron; 


PVR 


Pb 


four-momentum of outgoing 


I resonance. 



COMMON BFourVecs 

Contains the four-momenta in the rest frame of the resonance (RRF) where the positive 
^-axis is given by the boost direction j3. This COMMON block also contains the five four- 
momenta of the external particles. 

COMMON WvdWSpinors 

Contains the WvdW spinors of the external particles. 
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5.3 Subroutines 



A collection of some important subroutines: 

InitPhysCon Initializes the physical constants; 



ResProp Initializes the resonance properties; 

MC Generates an event using scheme I; 

MC2 Generates an event using scheme II; 

Boost Performs the boost to the RRF; 

StoreVar Stores the data of the generated event; 

Finish Calculates the cross sections after the generation loop 
and prints the results. 



5.4 Functions 

A collection of some important functions: 

RNF100 Generates a random number uniformly distributed in the 

interval [0, 1]; 

MatrixElement2 Calculates \M\ 2 using the expression in terms of invariants; 

SpinorMatrixElement2 Calculates the matrix element squared (and the density matrix) 

using the WvdW calculation with spinor inner products; 

FastSpinMat2 Calculates the matrix element squared (and the density matrix) 

using the WvdW calculation with traces of matrix products 
(preferred method for density matrix calculations); 

Norm Calculates the length of the spatial part of a four-vector. 



6 Comparison to other MC generators 

The results of the cross section calculations have been checked with the results of Galuga 



ljj . In this MC generator the same model is implemented in the Budnev || formalism. 
Within the errors we found complete agreement. Density matrix calculations could not 
be checked with any Monte Carlo generator, since only GaGaRes is able to calculate them 
at present. Only by extending the Galuga generator ourselves checks could be performed 
and gave agreement. 



7 Conclusions 

GaGaRes is a MC generator which is well-suited for the description of resonance produc- 
tion in two-photon physics. It offers several schemes and calculational methods and can 
be used for every required topology of the outgoing particles. The program is also able 
to generate density matrices which form an essential tool in the description of the decay 
of the resonances. 
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A Expressions for E \M\ 2 



In the expressions for the total matrix elements squared we have used the following 
ables to compactify the expressions 





= s n + s' n ±u n ±u 


0"2± 


= ss' ± uu', 


s n ± 


= s n ± s' n , 


U n ± 


= u n ± u' n , 


t n ± 


= t n ± t' n , 



with the convention that for n — 1 the n is not written. 

This leads us to the following expressions: 
1 5'o amplitude: 

\ e i-M(^o)i 2 = A{\ + «+) - («' + ^-) 2 - («' - ^ 2 -) 2 ; 

+m 2 [t + (s 2 + ui ) + 4<j 2 _£_ - 4tt'S + + 2s_w_t_ 
+16m^s + M + -64m^S + + 256m^} . 



3 P amplitude: 

\ E |-M( 3 A)| 2 = 4* {f (2<r 2 _ + tt'(s 2 + « 2 _ + 2*0) + Btf(-s + (a 2 - + tt>) 
+u + (a 2 „ - tt')) + 2t 2 t' 2 a 2+ + m 2 e [B 2 {t + Z 2 + - 4 ( x 2 _£_) 
+2tt'B(Y 2 _ - 2t+E+ + Att') - 4t¥ 2 (S+ - *+)] 
+m A e [AB 2 {Y?_ - 4t+S+) + 32tt't+B + lQt 2 t' 2 } +Q4m 6 e B 2 t + } . 

Here B is given by 

3 3 
B = h ■ k 2 + W 2 = -(s + s' + u + v!) + t + t' - 12m 2 = -E+ +t+- 12m 2 . 

3 P X amplitude: 

\ E |A^( 3 Pi)| 2 = {it [tt't 2 + - tet+E+ + \t 2+ ^2+ - 2tt% + 

+ (t 2+ - 4tt')<72+ + <7 2 -2 2 _ - 4tt's + M+ + t+(7 2+ S + 

+2t(sw' 2 + sV + s'u 2 + s' 2 u) + 2t'(sw 2 + s 2 u + sV 2 + s'V) 
+2t 2 (sw' + s'u) + 2t ,2 (sw + s'u')] - al_t 2+ 
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+m 2 [t 3 +S^ - 7tt'£+E 2 + - 26tt'(t(su' + s'u) + i'(su + s'u')) 
+8t 2+ (ss's + + uu'u + ) + 4(t 2+ + tt')(ss'u + + uu's + ) 
-22tt't + a 2+ - 12tt't 2 +E+ + 32t¥ 2 E+ - 2tt'E 3 + 
-6tt'( ss ' s + + uu'u + ) + 16t 2 t /2 t+ + 2tt'(s + u 2+ + u + s 2+ ) 
+4t 2 (su' 2 + s 2 u' + s'u 2 + s /2 w) + 4t' 2 (su 2 + s 2 u + sV 2 + s'V) 
-lOtt^sfa + siV + s't'u + s'tu')} 

+ 16mj [-t%s + u+ + 6tt't+E+ - t 3 +E+ - t 2+ (E 2 + + 4(j 2 +) 
+4tt't 2 - 2(t 2 (su' + s'u) + t' 2 (su + s'u')) 
+64m^ [2(i 2+ + + ^3+ - 5ft'*+] 

-256m^ + .} 



3 P 2 amplitude: 



|E |-M( 3 P 2 )| 2 = i [£ {$[-8E* - 32(72+2+ + 64 S+M+ E+ 

+192(ssV + uu's+)\ + #[-4E + + 128^+ - 16s+m+] 

+m 2 [-24E+] + ^|[-192S + - 256(x 2 + - 640s+«+] 

+#[-128E+] + m*[96] + #[3072E+] + #[512] + #[-8192] 

+^[E + - 2s+m+(s + + u + ) + 16cr|_ - 8s + u + a 2+ 

— 16(ss / -u 2 +m-u / s 2 +) — 4s 2 +m 2 + — A&ss'uu'] + [4E + — 8s+-u+E H 

-26o- 2 +E+ + 20(ss'w+ + uu's+)\ + [s+ + w+ + 8(x 2+ ]} 

+ + {#[2S 4 + - 128(T 2 _ - 32(x 2 _E 2 _] 

+^|[8E3 + 8a 2 _E_] + ml^T 2 ] + #[-32E^ - 128s+ M+ E+ 

+256(T 2 _E_] + #[-200E^ + 32s+w+] + m^[64E+] 

+ #[512E + + 1024s+m+] + #[1536E+] + m°[-256] 

+ #[-4096E+] + #[-4096] + $[8192] 
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+^[8<r 2 2 _E + ] + ^[-2<7 2 _]} + Sh24( S + u)(s> + «')] 
+^[-24(s + + «)] + i {#[8S 2 + + 64(s + + «) 
+ 128(sw' + s'u) + 32i'E + ] + #[-64E + - 48t'] + m 2 [48] 
+#[-640E+ - 256*'] + $[416] + #[2560] + ^[2E 3 + 
+2t'T? + - 8t'(s + u)(s' + v!) - 2A(ss'u + + uu's + ) - S(su'(s + u') 
+s'u(s' + u))] + ^[4E 2 - 4(s + u)(s' + u') - 12a 2+ ]} 
+i {#[8E^ + 64(s + «)(s' + u') + 128(sm + sV) + 32tE+] 
+#[-64E + -48t]+m 2 [48] 

+#[-640E+ - 256t] + $[416] + #[2560] + ^[2E 3 + 
+2£E^ -8t(s + w')(s' + w) 

-24(ss'w + + uu's+) - 8(su(s + «) + sV(s' + u'))] 
+ ^[4E 2 + - 4(s + + u) - I2a 2+ }} 

+W^2 {-^l[-56o-|_S + - 8cr 2 _E 3 _ - 8a 2 -(s + uu' - u + ss') 
+8a 2 _(s + u 2+ - u + s 2+ )] + #[-32<r 2 _ _ 8a 2 _E 2 _] + m 2 [16a 2 _E_] 
+#[416<r 2 _ + 160<r 2 _E 2 _ + 8E 2 + - 32s 2+ m 2+ ] + #[8E 3 + cr 2 _E_ 
— 32(s + -u 2+ + u + s 2+ ) — 64(ss'w + + uu's + )} + ml\8T? + + 64s + w + ] 
+ #[-128E 3 f - 512a 2 _E_ + 512s + «+E + ] + #[-64E 2 ] 
+m6[-384S + ] + #[512E 2 ] + m8[1536] + ^[2<j 2 _E 2+ + Aa 2 2 _a 2+ 
+4a 2 2 _s + u + ] + T ^[2 ( 7 2 _E + ] + [-4(7 2 _]} 

-128 S+M+ ] + #[-160E + ] + #[1536E + ] + #[640] + #[-4096] 
+ M^U) + h {#^'^] + #[24( S + u')(s' + «)] +#[-128t'E + ] 
+ #[512t']} + £ {#[8tE 2 + ] + #[24( S + «)(✓ + «')] + #[-128tE + ] 
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+$[512*]} + {$[64£ + + 32t+] + $[-144] + $[-512] + ^[2S 2 ] 
+^[12i: + + 6t + ] + [8]}' . 



1 D 2 amplitude: 

{^mQd^ = -j:\mQ s«)\\ (si) 

with 

c = °A (82) 

Cl 

and 

f "2r + 4^ J^-*+ 4W 2 J + ei~ W^) ■ (83) 

Note that for the 1 D 2 case the overall factor is included in the 1 S factor. 
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